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I. INTRODUCTION 

Recent advances in experimental investigations of 
single-electron molecular ionization in a low frequency 
strong laser field [T}^H] have created a demand for a 
theory of this phenomenon. There is a broad variety 
of theoretical approaches to this problem: molecular 
extensions of the analytical atomic strong-field meth- 
ods [TUH2U] . numerical methods explicitly incorporat- 
ing tunnelling (the Floquet approach [211 122] , the com- 
plex scaling method and the method of complex 
absorbing potentials [Mj)j numerical solutions of the 
time-dependent Schrodinger equation within the single- 
active-electron approximation |25H31j . numerical solu- 
tions of the time-dependent Schrodinger equation for 
two-electron systems and treatments based on 
the time-dependent density functional theory [38H4I] . 

As far as low frequency laser radiation is concerned, 
one can ignore the time-dependence of the laser and 
consider the corresponding quasistatic picture, which is 
obtained in the limit the laser frequency a; — > 0. In 
this limit, single electron molecular ionization is real- 
ized by quantum tunnelling. This approximation is valid 
from qualitative and quantitative points of view, and it 
tremendously simplifies the theoretical analysis of the 
problem at hand. Such single-active-electron approaches 
to molecular ionization, where an electron is assumed to 
interact with multiple centres that model the molecule 
and a static field that models the laser, are among the 
most popular. Analytical and semi-analytical versions of 
these methods, which are based on the quasiclassical ap- 
proximation [T4l - fT7| [T9] are indeed quite successful in 
interpreting and explaining available experimental data. 
However, these quasiclassical theories heavily rely on the 



assumption that the electron tunnels along a straight tra- 
jectory. The purpose of the current paper is to study the 
reliability of this hypothesis. 

Relying on the geometrical approach to many dimen- 
sional tunnelling by Hislop and Sigal which is 
a mathematically rigorous reformulation of the instan- 
ton method, we first introduce the notion of leading tun- 
nelling trajectories. Then, we analyze their shapes in the 
context of single-active-electron molecular tunnelling. It 
will be rigorously proven that the assumption of "almost" 
linearity of leading tunnelling trajectories is satisfied in 
almost all the situations of practical interest. Such re- 
sults justify the above mentioned models, and perhaps, 
open new ways of further development of quasiclassical 
approaches to molecular ionization. 

The rest of the paper is organized as follows: Section 
|ll]is a concise introduction to the Hislop and Sigal geo- 
metrical ideas and related topics. The proof of the results 
regarding the shapes of the leading tunnelling trajectories 
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are presented in Sec. HI We employ multiple spherically 
symmetric potential wells, which is the simplest type of 
molecular potentials, to estimate single electron molecu- 
lar tunnelling rates in Sec. |IV[ Leading tunnelling tra- 
jectories are numerically computed in Sec. |V] for model 
diatomic molecules of different geometries, and the rule 
of thumb on how to find the shapes of leading tunnelling 
trajectories is formulated. Conclusions are drawn in the 
last section. Finally, the Appendix contains the deriva- 
tion of the multidimensional generalization of the Landau 
method of calculating quasiclassical matrix elements. 



II. MATHEMATICAL BACKGROUND 

The instanton approach is one of the methods for the 
description of tunnelling [IS^ . It can be introduced as a 
result of application of the saddle point approximation 
to the modification of the Feynman integral obtained 
by performing the transformation of time t — >■ —ir to 
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"imaginary time" r (i.e., the Wick rotation). This tech- 
nique has turned out to be tremendously fruitful in many 
branches of physics and chemistry (see, e.g., Refs. [TTl - 

m)- 

We shall reiterate the main steps in deriving the in- 
stanton approach. Let us consider a quantum system 
with the Hamiltonian 



H 



-A/(2to) -}-C/(x), 



(1) 



where A is the n-dimensional Laplacian and x is an n- 
dimensional vector. The Feynman integral representa- 
tion of the propagator reads [5T] (atomic units are used 
throughout, unless stated otherwise) 



(x/le-^'olx, 



N / ^[x(i)]e'^[''(*)l. 



(2) 



5'fx] = / ^U,k)dt, J^(x,i) = U(x 

' 2m 



where the path integral sums up all the paths that obey 
boundary conditions x(0) = x,j and x(to) = x^-, and 
x(t) = dx.{t)/dt. After performing the Wick rotation, 
Eq. ([2| becomes 

{^f\e-"^° |x,) =N J f^[x(T)]e-^[^Ml, (3) 
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J_ /dxij) 
2m \ dr 



dr. 



where tq — Uq and S is called the Euclidian action. 
Hence, one can say that the transition from Eq. ^ to 
Eq. ([3]) is achieved by the following formal substitutions 

t^~iT, x{t)^x{T), x{t) ^ idx{T)/dT. (4) 

Comparing the actions S and S, one concludes that the 
motion in imaginary time is equivalent to the motion in 
the inverted potential. In other words, the actions S and 
S are connected by the substitution 



U ^ ~U, {E^^E). 



(5) 



The final step in the instanton approach is the applica- 
tion of the saddle point approximation to the Euclidian 
Feynman integral in Eq. ([3| assuming that tq — > 00. 

However, there is a long ongoing discussion jS^HSS] 
whether the instanton approach agrees with the quasi- 
classical approximation for tunnelling; some observations 
have been made that these two methods may disagree up 
to a pre-exponential factor. Furthermore, as it has been 
pointed out in Ref. |48j . the instanton approach in the 
formulation presented so far [substitutions Q] not only 
looks like a "highly dubious manoeuvre," but also gives 
no prescription for getting a correct pre-exponential fac- 
tor. In authors' opinion, such discrepancies mainly occur 
because the Feynman integral is just a heuristic construc- 
tion without a sound mathematical ground [5 6) . The ab- 
sence of strict rules of calculation of the Feynman integral 



makes impossible any definitive judgment of a particular 
result. Consequently, a natural question arises how the 
instanton method can be safely used and what the mean- 
ing of the substitutions Q and ([5| is. 

The mathematical physics community has reinter- 
preted the instanton approach rigorously (see, e.g., Refs. 
P^HISI [STHBH] and references therein), and the corre- 
sponding analysis answers both questions. Moreover, this 
rigorous interpretation is extremely useful because it can 
be implemented as an effective numerical method, which 
will lead to a clear physical picture applicable to a broad 
class of problems. We shall review briefly the above cited 
works since on the one hand, they are often unfamiliar to 
physicists, and on the other hand, they may be challeng- 
ing to read for non-specialists in mathematical physics. 

Historically, the first problem considered within such 
a framework was "how fast does a bound state decay at 
infinity?" [5Il [58] (see also Sec. 3 of Ref. [45]). Let 
us clearly pose the question. Consider the Hamiltonian 
([!]) as a self-adjoint operator on ^2iW^) - the space of 
square-integrable functions. A bound state wave func- 
tion ■0 S (K") is a normalizable eigenfunction of such 
a Hamiltonian, = Eip. Since the normalization inte- 
gral converges, the bound state wave function ip = ip{^) 
must vanish as ||x|| — > 00. Therefore, we want to deter- 
mine how this decay is affected by the potential U. This 
question can be answered very elegantly if we confine 
ourself to an upper bound on the rate of decay. 

To obtain this upper bound, we need to introduce first 
some geometrical notions. Let M be a real n-dimensional 
manifold (intuitively, M is some n-diniensional surface). 
The tangent space at a point x e Af, denoted by Tx{M), 
is a real linear vector space IR" that intuitively contains 
all the possible "directions" in which one can tangentially 
pass through x. A metric is an assignment of an inner 
(scalar) product to T^{M) for every x € M . 

Let X e K" and i,rj e T^{M). We define a (degener- 
ate) metric by 



(t»7>x = 2m([/(x)-i?)+(^,r,). 



(6) 



where {^,r]) = ^ ■ rj = S^iiji + . . . -I- ^nVn is the Euchdean 
inner product and /(x)+ = max{/(x),0}. Following the 
convention used in mathematical literature, we shall call 
metric Q as the Agmon metric. 

Having introduced the metric, we can equip the mani- 
fold M with many geometrical notions such as distance, 
angle, volume, etc. The length of a differentiable path 
7 : [0, 1] — ^ M" in the Agmon metric is defined by 



LAh) 



\i{t)h(t)dt 
2m [\u{'fit))-E]'/'\\jit)\\dt, (7) 



where ||^|| = ^/(^TO is the Euclidian (norm) length, and 
1 1 ^ 1 1 X = \/ a, ^)x- The path of a minimal length is called 
a geodesic. Finally, the Agmon distance between points 
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X, y € M", denoted by p£;(x, y), is the length of the short- 
est geodesic in the Agmon metric connecting x to y. 

Before going further, we would like to clarify the phys- 
ical meaning of the Agmon metric. Let us recall the Ja- 
cobi theorem from classical mechanics (see page 150 of 
Ref. 61J and page 247 of Ref. [SI]): The classical trajec- 
tories of the system with the potential f7(x) and a total 
energy E are geodesies in the Jacobi metric 



((1,77)), = 2m(ii;-C/(x)) + (1,77), 



(8) 



on the set {x e E"|C/(x) E} - the classical al- 
lowed region. The Agmon metric [Eq. ^] and the 
Jacobi metric [Eq. ([s])] are indeed connected through 
the substitution ([5]). By virtue of this analogy, we con- 
clude that the Agmon distance has to satisfy a time- 
independent Hamilton- Jacobi equation, also known as an 
eikonal equation. 



|VxPis(x,y)|2 = 2m(f/(x)-i?)H 



(9) 



where Vx/(x) = {df/dxi, . . . , df /dxn)- In fact, the Ag- 
mon distance is the Euclidean version of the reduced ac- 
tion [here, the adjective "Euclidian" means the same as 
in Eq. ([3])]. In other words, the Agmon distance is the 
action of an instanton. 

Now we are in position to recall upper bounds on a 
bound eigenstate of the Hamiltonian Q. First, under 
very mild assumptions on U (continuity, compactness of 
the classically allowed region, and absence of tunnelling, 
i.e., the spectrum of the Hamiltonian being only real), it 
has been proven [57] that for an arbitrary small e > 0, 
there exists a constant < Cg < oo, such that 



=2(1- 



0Pi=(x)|^(x)|2d"x< C, 



(10) 



where pe(x) 
means that ip{x.) 



P£;(x,0). Roughly speaking, result (10) 

■ (l-e)pB(x)- 



O (e"^^""^''^wj. However, this re- 
sult can be improved. For any small e > 0, there exists a 
constant < < oo, such that the following inequality 
is valid under additional conditions of regularity of the 
potential U 



|V'(x)| s^Qe-^i-^)''^^"). 



Analyzing Eq. (10) and Eq. (11), we conclude that 



the Agmon distance from the origin describes the expo- 
nential factor of the wave function. Further information 
can be found in Refs. [45l Ell [60] and references therein. 
We note that lower bounds on ground states can also be 
obtained by utilizing the Agmon approach |58j . 

We illustrate the power and utility of the upper bound 



(11) by deriving upper bounds for matrix elements and 



transition amplitudes in the Appendix. The former result 
is an estimate of the modulo square of the matrix element 
(i/'pl V \^|:q) , where ^pp and ipq are bound eigenstates of 
the Hamiltonian (|T| that correspond to eigenvalues Ep 
and Eq. It is demonstrated in the Appendix that for an 



arbitrary small e > 0, there exists a constant < Cg < cx3, 
such that 

\{i^p\V c, f l^2(x)e-2(i-^)[pB,(x)+pE,(x)]^„^^2) 



which could be interpreted as 

|2 



\{^p\vm\' 

= 0(^j V^2(x)e-2(l-^)[p«p(x)+PE,(x)]^n^^ (^3) 



Simplicity of the derivation of Eq. ( 13 ) does not imply 



its insignificance. On the contrary, Eq. (13) is a multidi 



mensional generalization of the Landau method of calcu- 
lating quasiclassical matrix elements [B^ (see also page 
185 of Ref. [Ml and Refs. [55] [SS]). To the best of 
authors' knowledge, such a generalization has not been 
reported before. To prove the one-dimensional version 
of the Landau method using analytical techniques (as it 
is usually done), one deals with the Stokes phenomenon 
(see, e.g., Ref. [67l); thus, the generalization to the mul- 
tidimensional case without too restrictive assumptions is 
not obvious. The Agmon upper bounds lead not only to 
quite a trivial derivation, but also to an intuitive physical 
and geometrical picture. 

Now we explain briefly how these geometrical ideas are 
generalized to the problem of tunnelling (interested read- 
ers should consult Refs. [i2M5] and references therein 
for details and further development). Let E be an en- 
ergy of a tunnelling particle. We denote the boundary 
of the classically forbidden region by 5_e. It is assumed 
that Se consists of two disjoint pieces and (i.e., 
Se = S];, U and 5^ n S"^ = 0) - the inside and 
outside turning surfaces, which are merely multidimen- 
sional analogs of turning points. Having introduced the 
concept of the Agmon distance, we naturally introduce 
two related notions: First, the Agmon distance from the 
surface S]^ to a point x, p_e(x, 5*^), as the minimal Ag- 
mon distance between the point x and an arbitrary point 
y & Se [more rigorously, p£;(x, S"^) = infygg^ p£;(x,y)]; 
second, the Agmon distance between the turning surfaces 
5^ and 5^, Pe{S^,S^), as the minimal Agmon dis- 
tance between arbitrary two points x e S'j^ and y G S"^ 



(11) [ pe{Se,S^) =inf^g5+pi5(x,S'£;)]. 



In a nutshell, and thus a bit abusing the formulation 
of the original result [H], we say that for an arbitrary 
small e > 0, there exists a constant c > 0, such that the 
tunnelling rate, F, (viz., the width of a resonance) in the 
quasiclassical limit [h — > 0) obeys 



rs$cexp[-2/3£(p£;-e)], 



(14) 



where < p^; < oo and I3ePe being the leading asymp- 
tote of pe{Sj^, S*^) when h ^ 0. However, the following 
interpretation of upper bound ( 14 ) is sufficient for our 
further applications: 



O e 



(15) 



4 



i.e., twice the Agmon distance between the turning sur- 
faces gives the leading exponential factor of the tun- 
nelling rate within the quasiclassical approximation. 

The Agmon distance between two points, pe{^, y), can 
be computed by solving numerically Eq. (|9| with the 
boundary condition 



PEiy,y) = 



(16) 



by means of the fast marching method [55H72) . More- 
over, having computed the solution, one can readily ex- 
tract the minimal geodesic from a given initial point x by 
back propagating along /9£;(x, y), where y is regarded as 
a fixed parameter; more explicitly, the minimal geodesic, 
g = g{t), is obtained as the solution of the following 
Cauchy problem [551 EI] 



k = -v^PEity), g(o) = x. 



(17) 



Such a geodesic can be interpreted as a "tunnelling tra- 
jectory." 

A brief remark on types of the solutions of Eq. ^ 
ought to be made. Generally speaking, an eikonal equa- 
tion admits a local solution under reasonable assump- 
tions, but a global solution is not possible in a general 
case owing to the possibility of development of caustics 
(see, e.g., Ref. ^TS')- Nonetheless, when we talk about 
a solution of Eq. ([9|, we actually refer to a viscosity 
solution because not only it is a global solution, but 
also it has the meaning of distance [701 IZI] which we 
originally assigned to the function pE- Another reason 
for employing only the viscosity solution of the eikonal 
equation is as follows: Writing the wave function as 
^'(x) — exp[— ^(/i; x)], the time-independent Schrodinger 
equation becomes 

\W^S{h;x)\^ - hAS{h;x) = 2m(C/(x) - E). 

Comparing this equation with Eq. ^ in the classical 
forbidden region, we conclude that 



Pe 



lim S{h; x), 



which is the definition of the viscosity solution of eikonal 
equation (|9| (see, e.g., page 540 of Ref. [75]). 

In fact, the fast marching method is an "upwind" finite 
difference method that efhciently computes the viscosity 
solution of an eikonal equation. Note, hence, that the 
fast marching method as well as the other ideas presented 
and developed in the current paper cannot be employed 
to study the influence of chaotic tunnelling trajectories 
(see Ref. |Hj and references therein) . Some implementa- 
tions of the fast marching method as well as the minimal 
geodesic tracing can be downloaded from Refs. [75H77] . 

The Agmon distance from the surface to a point, 
P£;(x, S*^), must satisfy Eq. Indeed, pe{x,S]^) is 

the solution of the boundary problem 



|VxPi=;(x,5^)|2 = 2™(C/(x)-^)_ 
PEiy,SE)=0, VyeS*^, 



(18) 



which can be solved by the fast marching method as well. 
Finally, the Agmon distance between the turning sur face s 
is computed as min^gg+ pe{^, S'^) after solving Eq. ( 18 1. 

The points b e 5^ and e e such that 



PEiSE,S;l^) = pE{h,e), 



(19) 



are of physical importance because they represent the 
points where the particle "begins" its motion under the 
barrier (b) and "emerges" from the barrier (e), corre- 
spondingly. Moreover, the minimal geodesic (17) that 
connects these points (g(0) = b and g(l) = e) is a tun- 
nelling trajectory which gives the largest tunnelling rates 
- the leading tunnelling trajectory. Note, however, that 
these points as well as the trajectories may not be unique 
in a general case. 

The idea of utilization of upper bounds to describe 
tunnelling is not new. Kapur and Peierls [^ [75] (see 
also Ref. |80J) have proposed as early as 1937 that even 
though many dimensional quasiclassical approximation 
is untractable in its original formulation, it still can be 
used to obtain the upper bound on the probability of 
transmission through a barrier. The geometrical ideas 
reviewed in the current section can be viewed upon as 
a reincarnation of the Kapur-Peierls approach with an 
important (and convenient for our applications) emphasis 
on the geometrical aspect of the method. 

It is also noteworthy that a power of the fast marching 
method in applications to tunnelling has already been 
recognized in chemistry within the context of the reaction 
path theory [8T1185j . Similarly to the current paper, the 
main object of interest of those studies is the reaction 
path, which is the leading tunnelling trajectory in our 
terminology. Nevertheless, the motivation for the usage 
of the fast marching method, presented in Refs. [8T1185j . 
is tremendously different from the geometrical point of 
view adopted here. 



III. MAIN RESULTS 

In this section, we shall follow a two step program. 
First, we consider tunnelling in multiple finite range po- 
tentials, where we prove that leading tunnelling trajec- 
tories are linear (Theorem [T]). Then, we reduce the case 
of multiple long range potentials to the previous one by 
employing the fact that a singular long range potential 
can be represented as a sum of a singular short range 
potential and a continuous long range tail [Eq. (36)]. 



Such a reduction allows us to prove that the leading tun- 
nelling trajectories are "almost" linear (Theorem [2]). We 
note that partitioning ( 36 ) was put forth by Perelomov, 
Popov, and Terent'ev [86-89 , and it is widely used for 
obtaining the Coulomb corrections in strong filed ioniza- 
tion (see Refs. [50^55] and references therein). 

Let us introduce some notations. Hereinafter, the di- 
mension of the space is assumed to be n ^ 2. The in- 
teraction of an electron with a static electric field of the 
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strength F is of the form Fxn {F > 0). dA denotes 
the boundary of the region A. The map, mirij;^ : M" D 
A — > M", selects a point x — min^^A e A that has the 
smaUest a;„ component among all the other points from 
A, assuming that A has such a unique point. The pro- 
jection Px of the point x = (xi, X2, ■ ■ ■ , Xn) is defined as 
Px^{xi,...,x„^i,E/F). 

Theorem 1. We study single electron tunnelling (—oo < 
E < 0, F > 0) in the potential 



K 



C/(x)-^V,(||x-R,||) + Fx, 



(20) 



Let us assume that 

1. Vj : {0,Rj) (-00,0) and Vj : {Rj, +oo) -> {0}, 
i?j > 0, J = 1, . . . , K , are differentiate on (0, Rj) 
and strictly increasing functions, such that Vj{Q) = 
— oo and Vj may have a jump discontinuity at the 
point Rj . 

2. supp V,- = {x e K" I V,(||x - Rjll) 7^ 0} is the sup- 



port of the potential Vj(||x 



R,7 



suppVfc n 



supp Vj = 



Vfc 7^ i and supp Vj n 



{XG 



c„ ^ E/F} = 0, J = 1, 



Introduce (\j = minj,^ 9supp Vj , Pj = 

va.inx„S]^{j), S^{j) is defined in Eq. (22) 
If there exists N , such that 



(21) 



Then, the leading tunnelling trajectory is unique and lin- 
ear, and it starts at the point pjy and ends at PRn , 

Proof. The boundary of the classically forbidden region 
is defined by the equation J7(x) = E. Consider two cases: 

First, if X ^ UjLiSuppVj then according to assump- 
tion 2, the equation C/(x) E simply reads Fxn — E, 
and thus its solution defines the outer turning surface 
5+ = {x e M" I x„ = i?/P}. One can see now that the 
projector operator P projects a point onto S*^. 

Second, if x G supp Vj and Vj is continuous at the point 
then the equation reads V,(||x — R^ ||) -I- Fxn = E. 
To proof that the set 

Se{3) - {x G suppy, I y,(||x - R,||) + Fx^ = E} (22) 

is not empty, we construct the function /j(x) = Vj(||x — 
RjII) -I- Fxn — E. Since fj^Rj) = — oo, we can find a set 
Aj C supp Vj located close to Rj , such that fj (x) < for 
all X G Aj-; correspondingly, since according to assump- 
tion 2, Xn > E/F, there exists the set Bj C suppV, 
of points close to the boundary of supp Vj for which fj 
is positive. In fact, Aj and Bj can be constructed such 
< ||y - RJ, Vx G Aj and Vy G Bj. 



that 



Rn 



Therefore, the intermediate value theorem guarantees 



that S^{j) 7^ and it "lies between" Aj and Bj. Fur- 
thermore, the inner turning surface is 5^ = UjLi ^eU)^ 
and S'^(j) n S'^(fc) = 0, Vj ^ k. (Note that the strict 
monotonicity of Vj{x) assures that the set 5^(j) is con- 
nected.) Whence, 



Pe{Se,S^) 



mm 

j 



{pe{Se{j),S+)} 



(23) 



Eq. (23) means the reduction of the many centre case 



to the singe centre case under the assumptions made. 
Needles to mention that such a reduction tremendously 
simplifies the analysis. 

The same conclusions are valid if the jump of the func- 
tion Vj at Rj is not too large, so that the equation 
Vj (||x — Rjll) -f Fxn = E has solutions for x G supp Vj-. 
However, if the jump is too large, i.e., this equation does 
not have solutions from the support of the potential, then 
it is natural to set S'^(j) — cJsuppl/,-. 

Consider the single centre case - single electron tun- 
nelling in the potential Uj{x) = V,(||x — Rj\\) -\- Fxn- 
We shall show that this potential is axially symmetric. If 
X = {xi, . . . , a;„), then we introduce IIx = (xi, . . . , Xn-i)- 
We can then symbolically write x — (IIx, a:„). Using this 
new notation, we obtain 



Hx-nRjl 



\R 



JinJ 



Fx^^A) 



where [a]^ denotes the n**^ component of the vector a. It 



is readily seen from Eq. ( 24 ) that the potential Uj (x) is 



invariant under transformations that do not change Xn 
and arbitrary (n — 1) dimensional (proper and improper) 
rotations of the vector IIx about the point HRj. The 
only invariant subspace of K" under such transformations 
is the line {(IIRj, a;„) | a;„ G M}. 

Since both regions S'^(j) and are shape invariant 
under the axial symmetry transformations, we may ex- 
pect that the shortest geodesic connecting these regions 
ought to be shape invariant as well. Thus, one readily 
concludes that the leading tunnelling trajectory should 
be linear and should connect the points pj and PRj 



PEiPj,PR,)^PE{SEU),S+) 



(25) 



since no other geodesic that connects S'^(j) and S*^ is 
shape invariant with respect to the axial symmetry trans- 
formations. Below we shall present a formal version of 
this derivation. 

Foremost, we demonstrate that the operation miria;^ 
is defined on the set S~{j), viz., that there is a unique 
point of S~ {j) that has the smallest component x„. Em- 
ploying the method of Lagrange multipliers and taking 
into account the symmetry of the potential, we construct 
the function 



+ X Vj [\]c^+{xn 



\R 



JinJ 



Fx„ - E 



.(26) 
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The condition d^i/dc = leads to c = 0. Therefore, 
Pj — min.j;^ S]^{j) = (nRj, y), where y being the mini- 
mal solution of the equation 



vA\y-i^j]n\)+Py = E. 



(27) 



Moreover, Ppj = PRj = Pqj. 

Eq. (27) must have two distinct solutions yi_2 (j/i < 
2/2)- 2/1 (^2) corresponds to the point from S]^{j) with 
the minimum (maximum) Xn- Additionally, since E — 
Fy,>E- Fy^ ^ V, {\y, - [R,] J) > V, (|y2 - [R,]J), 
we obtain 



Vj = \vi - [R-j]„| > I2/2- [Rj]„| 



(28) 



To find the maximum of the function |jx — Rj|j on the 
set S^{j) within the Lagrange multipliers method, we 
introduce the function 



^2(a;„,c, A) = \fcA+Jx 
+A 



fR, 



fR,- 



Fx..r - E 



and the fact that 
lat the maximum of 



Taking into account inequality (28 
8^2/ dc = =4> c = 0, we conclude t 
the function ||x — Rj|| on S^{j) is reached at the point 
p,. 

Let Sj{z) denote a sphere of the radius z centred at 
Rj, Sj{z) = {x G E" I llx - Rjll = z). Consider a se- 
quence of spheres {Sj {rjj + k[Rj — ?7j]/W^)}J^o , where 
Sj{Rj) ~ i9suppVj and r]j was introduced in Eq. (28 1. 
Now pick a sequence of points, {7(fc/VF)}^g, such tnat, 
■j{k/W) £ Sj {r]j + k[Rj - r]j]/W), k ^ 0,...,W. We 
assume that this sequence is a discretization of some dif- 
ferentiable path 7 : [0, 1] ^ M". According to Eq. 
the sums. 



w 



2mY,Ju,h{k/W})~E 



fe=0 



x||7([fc-|-l]/W^)-7(fc/M/)||, (30) 



where we set 7(1 + l/W) = 7(1), obeys the property 
limiy_).oo ^wil) = ^Ail)- Introduce a path: 



g{t) = Pj+t [qj - Pj] . 



(31) 



Since Vfc, g{k/W) e Sj {ri, + k[Rj ~ iy]/W), 
hik/W)]^ > [g{k/W)l and V,i\\gik/W) - R,||) = 
VMk/W) - R,||) ^ U,{^ik/W)) ^ UMk/W)). 
Moreover, WfHk + 1]/W) ~ j{k/W)\\ ^ \\gi[k + l]/W)- 
g(fc/W^)||. Therefore, 



(32) 



9supp Vj and must be a straight line that starts at 
Qj and ends at Pqj because the potential between these 
two regions is merely ^(x) = Fxn. 

To finalize the proof, we shall "backward propagate" 
the leading tunnelling trajectory starting from the outer 
turning surface S^. Let p(x, y) denote the Agmon dis- 
tance between two points for the potential V^(x) — Fxn- 
Then, it is easy to demonstrate that 



Pe{x, Px) (2/3)V2toP||x - Pxf/^. 



(33) 



The plane T{c) = {x G M" | a;„ = c} is a sur- 
face of constant Agmon distance [Eq. (33)], such that 
Pe{Te/f,S^) = and pe{T{c) 



creasing function of c. Since Hq^y 
PRatII - IIpjv - qwll < ||q, - PR 



is a strictly in- 
PRatII = IIpat - 
J - ^ Vj ^ N, condi- 

tion ( [2l| ) guarantees that increasing c the plane T(c) will 
"hit" the boundary of suppV^r at the point qAr. (Note 
that pE {T{c),S+) = pE {T{c),S+), E/F<c< [qw]„.) 
Moreover, the following follows from Eq. (21) 



•(29) {x G M" I [qAr]„ sc: x„ < [pat] J n supp V, = 0, Vj ^ N, 



which means that iV"^ centre is isolated from all the 
other. Therefore, the shortest geodesic must connect the 
point qAT to the point pat. □ 

Corollary 1. Consider a single electron tunnelling in the 
potential \2(^ , such that assumptions 1 and 2 of Theorem 
[7] are satisfied, then the leading trajectory is linear (but 
may not he unique). 

Proof. This corollary follows from the straightforward 
generalization of the idea of backward propagation. □ 

Theorem 2. We shall study single electron tunnelling 
(—00 < E < 0, F > 0) in the potential 



K 



C/(x) = ^r,(||x-R,||) + Pa;„. 
that 



(34) 



1. 



(0, +00) — !■ (—00, 0) are differentiable on 
(0, -l-oo) and strictly increasing functions, such that 
YjiO) = -00 and Yji+oo) = 0. 

2. The boundary of the classically forbidden region 
consists of two disjoint pieces - the internal turning 
surface S"^ and the outer one S^. Furthermore, 

Se = [Jf=iSEU), s^U) n s^ik) = 0, Vj ^ k, 

where each S'^(j) encircles Rj^. 



Since J^wh) = ^wis) ^ l{k/W) = g{k/W), k = 
0, . . . ,W — 1, WW, we conclude that path ( [3l] ) is indeed 
the shortest geodesic that connects the regions S'^{j) and 
3supp Vj . By the same token, the geodesic connecting 



^ The verb "encircle" should be understood in the following sense: 
A piece of the inner turning surface, S^{j) = dCA{j), is a 
boundary of the classically allowed region, CA{j), associated 
with centre j, such that Rj £ C'A{j). 
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3. B{j) n B{k) = 0. Vj 7^ fc, and B{j) n 5*+ = 0, 
Vj, w/iere B{j) = {x e M" | ||x - Rj|| r^} 6e- 
mi? the ball of radius rj centered at Rj . Here 
Tj = max {||x — Rjll I X € is the "radius" 

Then, the leading tunnelling trajectory (may not be 
unique) is linear up to a term of 0{\) as X ^ Q, where 
A = maxj {|-^(Aj)|} and 

A, = min (^-j + l mm^. " Rfell " ^k} , d,^ . (35) 

Here, dj — min | ||x — R^ || | x G S'_^} is the Euclidean 
distance from Rj to . 

Proof. We introduce two auxiliary functions 








< a; < A 
X ^ Aj, 

Os^x < A 
X ^ Aj. 



J' 



J' 



One evidently notices that 



(36) 



where V^ji\x) is a singular short range potential and 

Vi^\x) being a long range tail. The purpose of such 

a partition is to make V^l\x) satisfy assumption 1 of 

Theorem jlj and produce Vi^J\x) that obeys the following 
upper bound: 



\v,f{x)\^\y,{A,)\^x, 



We analyze the length of a curve in the Agmon metric 
[Eq. Since 



VU{X} - E : 



K 



\ Y.v}^\\\^~'R,\\)+Fxn-E + OiX) 



K 



E ^shiW^ - II) + F^n -E + 0{X), 



under the assumption that A ^ 0, we have reduced the 
initial situation to the case of single electron tunnelling 
in the potential 



K 



C/,,(x)=^Fi^")(||x-R,||) + ^^x, 



(37) 



The parameter Vj can be calculated by means of the method of 
Lagrange mu ltipliers as it was shown in the proof of Theorem [T] 
[see Eq. E9\]. 



Let us now utilize assumption 3 to show that 



(38) 



> rj ; on the 



IIR,- — Rt 



-rk > 



Indeed, on the one hand, B{j)r\S]. 
other hand, B{j)r]B{k) ^ 0, Vj ^ k, 

Furthermore, we shall demonstrate that the definition 
of Aj [Eq. (35l] assures that assumption 2 of Theorem [l] 



for the functions V^l^ {x) holds. According to Eq. (35), 



A, s^dlR, -Rfc||-rfc + r,)/2, j ^ fc; 



|R-j— Rfcl 



supp v^l^ nsupp yj^-' 



suppyj;^^ n 



hence, Aj+Aj. < 

From Eq. (35), we also obtain Aj ^ dj 

= 0; thus, the outer turning surface for the potential 
Ush{y^) should be {x e M" | Xn = E/F}. 

Finally, we have proven the theorem because the po- 
tential Ushi^) satisfies all the assumptions of Corollary 
[D □ 

Physical clarifications of Theorems [l] and [2] are due. 
Assumption 1 of Theorem[T]physically implies that Vj are 
attractive, singular, spherically symmetric short range 
potentials. Assumption 2 of the same theorem requires 
that the potentials do not merge, i.e., their ranges do 
not overlap. This condition connotes that the classi- 
cally allowed regions associated with the centres Rj [their 
boundaries are S]^{j)] do not overlap as well. The lat- 
ter statement is proven in Theorem [T] The statement 
of Corollary [T] can be rephrased as follows: leading tun- 
nelling trajectories for a system of non-overlapping, at- 
tractive, singular, short range potentials are linear. How- 
ever, if the additional condition (21) is satisfied then The- 
orem [T] not only guarantees the uniqueness of the leading 
tunnelling trajectory, but also provides the initial and 
final points of the trajectory. Assumption 1 of Theo- 
rem [2] means that yj are attractive, singular, spherically 
symmetric long range potentials that vanish at infinity. 
Assumptions 2 and 3 of Theorem [2] require the same non- 
overlapping condition for the classically allowed internal 
regions mentioned above. Physically, Theorem [2] says 
that leading tunnelling trajectories for a system of sev- 
eral such potentials are "almost" linear, and a deviation 
from being strictly linear is caused by vanishing long tails 
of the potentials; thus, the larger the distance between 
the centres, the smaller the deviation. 

In a nutshell, all these results have been achieved be- 
cause the multi centre (i.e., molecular) potential is repre- 
sented as a sum of spherically symmetric potentials, and 
such conclusions regarding the shape of the trajectories in 
the single centre (i.e., atomic) case are quite expectable 
owing to the axial symmetry. 

An important case not covered by the theorems is 
the case of overlapping potentials that physically corre- 
sponds to valence electrons, which form chemical bonds 
and have a low ionization potential and are delocalized 
over a molecule. This case as well as the issue of unique- 
ness of the trajectories will be scrutinized in Sec. |Vj 
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IV. THE APPLICATION OF SPHERICALLY 
SYMMETRIC POTENTIAL WELLS TO SINGLE 
ELECTRON MOLECULAR TUNNELLING 

The simplest type of model molecular potentials that 
allows for full analytical treatment is of type ( 20 ) where 



— oo 




< a; < r 

X > Ti 



J' 



(39) 



It is assumed that S^{j) — 9supp = {x e M" | ||x|| = 
rj}. (Strictly speaking, these potentials are not governed 
by Theorem [I]) Evidently, the leading tunnelling trajec- 
tories are linear, and moreover, the following equality is 
valid 



Pe{Se,S%) ^\-nin{pE (qj,Pqj)}, 



(40) 



where — ra.\nx,^S^{i) and pE was defined in E g. (|33| . 
Let us estimate the tunnelling rates within Eq. (|l5| for 
the two dimensional system of two equivalent centres of 
type (39) (see Fig. [T]). A straightforward geometrical 
derivation, using Eqs. (IT5|, (33 1, and (40 1, shows that 



^2 




Se(2) 




-I' 


/ 

\ / / 






r + R/2 







E/F 





FIG. 1: The geometry of a two centre model employed to 
obtain Eq. (41 1. Grey colour denotes the classically allowed 
regions. 



Fcxexp^-— [Fi?(I 



I cos 6*1) - 2E\ 



3/2 



(41) 



where R is the distance between the potential wells (i.e., 
the bond length of a model molecule) and Q is the angle 
between the field and the molecular axis. The obtained 
angular dependent rates are plotted in Fig. [2] 

According to Eq. (15), rates obtained within the geo- 



metrical approach do not account for an initial molecular 
orbital. This technique provides solely the contribution 
from the shape of the barrier, hence, the name - the "ge- 
ometrical approach." An advantage of such a method is 
that it reduces the calculation of tunnelling rates to a 
rather simple geometrical exercise. 



180 




270° 



V. NUMERICAL ILLUSTRATIONS 

In order to illustrate the results of Theorem [2] and also 
draw some conclusions beyond Theorem [2] we shall cal- 
culate the shapes of leading tunnelling trajectories for 
different situations within the numerical scheme sketched 
in Sec. |ll] To achieve a good accuracy, we compute the 
viscosity solutions of eikonal equation ([9]) by means of 
the second order multi-stencil fast marching method [72] . 
The following two-dimensional model potential is used for 
a diatomic molecule in this section: 



2 



+FX2 



-1/2 



(42) 



FIG. 2: The polar plot of the normalized tunnelling rates for 
the two centre model [Eq. (41 1] as a function of the angle 9. 
Chosen values of the parameters are F = 0.01 (a.u.), R = 2 
(a.u.), and E — —0.5 (a.u.). 



where the first atom is centred at {^[^\^2^^) and the sec- 
ond atom is at (^p-* , ^2^-* ) , 

^ _(i?/2)sin6i, ^2^' =4 - (i?/2)cos6l, 
^ = (R/2) sin 0, ^2^' = 4 + (i?/2) cos 6, 
= 0.05 (a.u.), = -0.5 (a.u.), m = 1 (a.u.).(43) 

Regarding the definitions of the interatomic distance R 
and the angle 9 between the molecular axis and the ex- 
ternal field F, see Fig. [l] 
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Solutions of the eikonal equation and leading tun- 
nelling trajectories for identical non-overlapping long 
range potentials ( 42 ) are pictured in Fig. [S] Figure [Sj^a) 
illustrates the non-uniqueness of leading tunnelling tra- 
jectories in the case of non-overlapping potentials. The 
reason of this non-uniqueness is the mirror symmetry of 
the potential with respect to the axis xi = 0, as a result, 
the probabilities of tunnelling along each trajectory coin- 
cide. Note that the leading tunnelling trajectories in Fig. 
3 ^a) are almost linear. Results shown in Figs, [sjjb) and 
3 c) are in full agreement with the statement of Theorem 

2 

Leading tunnelling trajectories for identical overlap- 
ping long range potentials are shown in Fig. |4] The 
shapes of all these trajectories are almost linear as well. 
Non-uniqueness of leading tunnelling trajectories in the 
case of overlapping potentials is demonstrated in Fig. 
|4|^a). However, if the interatomic distance is further de- 
creased, one observes in Fig. |4]jd) that the previous two 
distinct trajectories merge into one restoring the unique- 
ness of the leading tunnelling trajectory. 

The shapes of leading tunnelling trajectories in the 
case of non-identical long range potentials are shown in 
Fig. |5] The interatomic distance is chosen such that 
one observes a transition between non-overlapping and 
overlapping cases by simply rotating the molecule. We 
clearly see in Fig. [5] that leading tunnelling trajectories 
are almost linear. 

The results presented in all these figures can be 
summed up in the following rule of thumb on how to 
find the leading tunnelling trajectory: The final point (or 
multiple points when more than one leading tunnelling 
trajectory is possible) is near the point on the outside 
turning curve (surface in the n ^ 3 dimensional case), 
, with the smallest value of X2 {Xn in the n dimensional 
case). The initial point is usually near the point on the 
inside turning curve, S'^, with the smallest value of X2- 
The tangent vector to the leading tunnelling trajectory 
at the initial and final points tend to be perpendicular to 
the inside and outside turning curves, respectively. 

Note that the rule of thumb is valid for polyatomic 
molecules and for an arbitrary number of dimensions. 



VI. CONCLUSIONS AND DISCUSSIONS 

Having introduced the leading tunnelling trajectory as 
an instanton path that gives the highest tunnelling prob- 
ability, we have proven that leading tunnelling trajec- 
tories for multi-centre short range potentials are linear 
(Theorem [T]) and "almost" linear for multi-centre long 
range potentials (Theorem [2] and the rule of thumb from 
Sec. 0. 

The fact that the leading trajectories for long range po- 
tentials are not straight lines is of vital importance. As 
in the atomic case EH-ES , this deviation is cru- 





-8 -6 -4 -2 2 4 6 



xi (a.u.) 




FIG. 3: (Color online) Viscosity solutions of eikonal equation 
(|9| and leading tunnelling trajectories in the case of two iden- 
tical non-overlapping long range potentials, which are given 
by Eq. (liSl) with Zi = Z2 = I and R = 6 a.u. White colour 



denotes the classically allowed regions. Dashed white lines are 
leading tunnelling trajectories. The solutions of the eikonal 
equation is represented by linear scale colour ramps from 
black (minimum) to bright color (maximum), (a) 6 — 90°; 
(b) 6 = 45°; (c) e = 0°. 
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-4-2 2 

(a.u.) 



FIG. 4: (Color online) Viscosity solutions of eiltonal equati on ((ot and leading tunnelling trajectories in the case of two identical 
overlapping long range potentials, which are given by Eq. ( |42[ ) with Z\ = Z2 = 1. White colour denotes the classically allowed 
regions. Dashed white lines are leading tunnelling trajectories. The solutions of the eikonal equation is represented by linear 
scale colour ramps from black (minimum) to bright color (maximum), (a) Q — 90° and R — 5 a.u.; (b) 9 = 45° and R — 5 a.u.; 
(c) = 0° and = 5 a.u.; (d) 9 = 90° and R = 3 a.u. 



cial for a quantitative treatment [TTJ [MHTT] , and some- 
times even for a qualitative analysis, because it leads to 
the correct pre-exponential factor of ionization rates that 
describes the influence of the Coulomb field of nuclei. 
However, Theorem [2] as well as illustrations presented in 
Sec. |V] suggests that the exact shape of the trajectory 
can be obtained as a perturbation of the trajectory in 
the modified potential where the Coulomb potential is 
substituted by the finite range one. This is a part of the 
celebrated Perelomov-Popov-Terent'ev (PPT) approach 
widely employed in the literature for analytical 
calculations of atomic Coulomb corrections. Neverthe- 
less, the PPT method requires matching the quasiclassi- 
cal wave function of an electron in the continuum with 
the bound (field-free) atomic wave function. This step is 
a stumbling block for generalization of the PPT approach 



to the molecular case (for the suggestion of a solution to 
such a problem see Refs. [ini US]). Theorem [2] in fact 
offers a solution to the problem of matching. According 
to Theorem [2] matching should be done on spherical sur- 



faces of radii [Eq. (35)] centred at the nuclei. This is 



an alternative technique to the method developed in Refs. 
[ini HH] ; however, it is applicable only to well separated 
core electrons and may not be applicable to delocalized 
valence electrons. 

It has been suggested in Ref. [97] that molecular pho- 
toionization in the tunnelling limit may act as a scanning 
tunnelling microscope (STM). Since rotating a molecule 
with respect to a field direction is analogous to moving 
the tip of an STM, then the observed angular-dependent 
ionization probability should provide information for a 
molecule similar to the position dependence of the tun- 
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(b) 



1 

J 



-8 -6 -4 -2 2 4 6 8 

xi (a.u.) 



(c) 




FIG. 5: (Color online) Viscosity solutions of eikonal equation 
Q and leading tunnelling trajectories in the case of two non- 
identical long range potentials, which are given by Eq. (42 1 
with Z\ = 0.4, Z2 = 1, and R — 3.9 a.u. White colour denotes 
the classically allowed regions. Dashed white lines are leading 
tunnelling trajectories. The solutions of the eikonal equation 
is represented by linear scale colour ramps from black (mini- 
mum) to bright color (maximum), (a) 8 = 90°; (b) 6 — 45°; 
(c) e = 0°. 



nelling current in the STM. We point out that there 
is a resemblance between such a descriptive comparison 
and our results. As we have demonstrated (the rule of 
thumb and the backward propagation of the leading tra- 
jectory), the leading tunnelling trajectory starts at the 
atomic centre that is the closest to the barrier exit (i.e., 
the outer turning surface); hence, the qualitative similar- 
ity of molecular tunnelling with the STM. Nevertheless, 
one must bear in mind that an electron cannot tunnel 
along a path because Heisenberg's uncertainty principle 
requires a wave packet with a finite lateral extension. It 
is demonstrated in Ref. [98j how the tunnelling current 
is drawn out of a bound state along the direction of the 
external field, spreading out also somewhat in the or- 
thogonal direction, and giving rise to a tunnelling "spot" 
of approximately a Gaussian shape. Such wave packets 
are important for the interpretation of experiments (see, 
e.g., Refs. [UlIlT]). 

The demonstrated simplicity of the shapes of leading 
tunnelling trajectories, in fact, may encourage future de- 
velopments of analytical models of molecular ionization. 
Nevertheless, the geometrical approach has a fundamen- 
tal limitation - it does not account for effects of molec- 
ular orbitals, and there is no a priori way of including 
these effects. In spite of that, one may always attempt 
to introduce such corrections in a heuristic manner, e.g., 
multiplying the geometrical rates by a Dyson orbital. 

In the current paper, we modelled a molecule by a 
single-electron multi-centre potential, hence discarding 
effects of electron-electron correlations. However, the ge- 
ometrical approach to tunnelling reviewed in Sec. [TT]can 
account for these effects after an appropriate adaptation 
presented in Ref. [33]. Intuitively speaking, according 
to such a method, the leading tunneling trajectory of the 
system is selected such that the minimum number of elec- 
trons are displaced during tunnelling. More importantly, 
the fast marching method still can be utilized to obtain 
this leading tunneling trajectory. Since correlation dy- 
namics of electrons plays an important role in molecular 
ionization leading to interesting novel effects [HI HSl 1551 - 
157] . applications of the geometrical ideas developed in 
Ref. [43] to molecular ionization should be the aim of 
subsequent publications. 
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Appendix A: Upper bounds for matrix elements and 
transition amplitudes 

The current section contains simple derivations of the 
muhi-dimensional generahzation of the Landau rule for 
calculation of the quasi-classical matrix elements between 



bound states [Eq. (A4)] as well as estimates of perturba- 



tion theory transition amplitudes [Eqs. (A6) and (A9)]. 
Note, nevertheless, that we do not employ these results in 
the paper. The purpose of these derivations is to demon- 
strate methodologically the utility of the Agmon upper 
bounds for bound states reviewed in Sec. |TT] as a pre- 
lude to the Agmon geometrical ideas used to describe 
tunnelling. 

For the sake of simplicity, the argument x will be 
omitted in some equations below. Throughout this Ap- 
pendix, we assume that the Agmon upper bounds [57] 
for bound states (-0„) are valid, i.e., Ve > 3c„ = c„(e), 
< c„ < oo, such that 



(Al) 



where pn = Pe„ ■ 

Let us choose an arbitrary e > 0. Employing the 
Schwartz inequality and assumption (All, we obtain 



= (l-e)(Pp-(-pg) 



Blg4ie')4{e') / e-2(.-.')(Pp+P,)rfx, (A2) 



where e > e' > and 



|y|ie-i(i-E)(Pp+p,)rfx. 



(A3) 



The integral / exp[—2{e — e'){pp + Pq)]dx converges for all 
p and q. Therefore, we have proven that Ve > 3c = c(e), 
< c < cx), such that 



tp*Vtp„dx 



(A4) 



which is the same as Eq. ( 12 ) 



Now let us study the problem of estimating of transi- 
tion amplitudes defined by means of the time dependent 
perturbation theory. Hereinafter, we assume that a quan- 
tum system under scrutiny has no continuum spectrum, 
and we shall manipulate with all the series and integrals 
heuristically - assuming that they all converge, or alter- 
natively, assuming that they are over a finite range. We 



illustrate our idea by estimating the second order ampli- 
tude since generalization to higher orders is evident. 

The second order transition amplitude within the time 
dependent perturbation theory reads 

= - f' y"*'dt'y"dxdxv;,„(x')e-*^^-(*^-*') 

xF(x')i^(xY|xt)l/(x)i/^„,(x)e-''^-(*-*'), (A5) 

where all the V^'s are eigenstates of the system and K is 
the propagator, which can be written as 



if(xY|xi)-^^„(x')C(x)e-* 



-iE„{t'-t). 



whence, 



|i^(xY|xt)K^|^„(x')^„(x)|. 



Using such a simple estimate as well as inequality ( Al I, 
we obtain 



(A6) 



where M = CinCfin max„ {c^} /2, < M < oo. 

However, there is no need to confine ourself to the case 
when the initial and final states are eigenstates. The 
same idea applies to the general case of the initial ((/>m) 
and final {4'fin) states being represented as linear expan- 
sions in the basis of the bound eigenstates, 

4>in = '^{'^n \<t>in) V'n, <l>fin = '^{'^n |0/m) "iAn- (A7) 
n n 

Let us found an upper bound for the first order transition 
amplitude, which is as follows 

= -I r dt J dx^(^;,„|Vn)C(x)e-'^"(*^-*) 

xT/(x)(V„' |0.„) ^„-(x)e-^^"' (A8) 
Whence, we readily obtain 



^ CnC„' |?A„) \(f>in) \B'i^^, 



n,n 



(A9) 



where M = max„^„' {c„c„' | ((/)/i„ IVn) (V'n' I'/'m) |}, < 
M < oo. 
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